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Abstract 

The  purpose  of  this  paper  is  to 
compare  two  distinct  methods  to 
extract  a  linear  state-space'  model 
about  a  reference  flight  condition 
from  a  nonlinear  simulation.  The 
frequency  response  method  uses  a 
time  history  input  which  contains 
the  frequencies  of  interest  to  drive 
the  simulation.  The  frequency  input 
and  the  output  of  the  simulation  are 
transformed  to  the  frequency  domain, 
and  the  desired  frequency  responses 
of  the  simulation  are  calculated.  A 
linear  model  is  then  fit  to  the 
frequency  responses  using  system 
identification  techniques.  The 
perturbation  method  extracts  a 
linear  model  by  perturbing  the  model 
states  and  inputs  from  the  reference 
flight  condition  and  calculating  the 
resulting  model  coefficients.  Both 
methods  were  used  to  extract  a 
fourth  order  longitudinal  state- 
space  model  from  the  V-22  full 
nonlinear  simulation.  The  time 
history  responses  and  system 
matrices  of  the  extracted  models 
were  compared.  The  comparison  showed 
that  both  methods  are  effective 
means  to  reduce  a  nonlinear 
simulation  to  a  linear  state-space 
model . 


Nomenclature 

ax  X  body  axis  acceleration 

( ft/sec^ ) 

az  z  body  axis  acceleration 

( f t/sec^ ) 

A  State-space  system  matrix 

B  State-space  control  matrix 

C  State-space  output  matrix 

Cxy  Coherence  function 
D  State-space  output  control 

matrix 

Chirp-z  transform  starting 
frequency  (rad/sec) 


^2  Chirp-z  transform  ending 

frequency  (rad/ sec) 

H  Transfer  function 

Ux  Length  of  windows 

bssw  Length  of  sum  of  sine  waves 

signal 

Nw  Number  of  windows 
Pxx  Input  power  spectral  density 
(dB) 

Pxy  Input -output  cross  spectral 
density  (dB) 

Pyy  Output  power  spectral  density 
(dB) 

qb  pitch  rate  (deg/sec) 
t  Time  (seconds) 

u  State-space  input  vector 

ub  X  body  cucis  velocity  (ft/sec) 

wb  z  body  axis  velocity  (ft/sec) 

X  State-space  state  vector 

X  Input  signal  Fourier 

Coefficients 

y  State-space  output  vector 

Y  Output  signal  Fourier 

Coefficients 

5e  Elevator  position  (deg) 

5ion  Longitudinal  stic)c  position 
(in) 

At  Simulation  time  step  (sec) 
Aco.-zt  Chirp-z  transform  frequency 
resolution  (rad/sec) 

AWssw  Sum  of  sine  waves  signal 

frequency  resolution  (rad/ sec) 
A,  eigenvalue 

A  dynamic  mode 

(0  Frequency  (rad/s_j) 

0  Pitch  attitude  (deg) 

*  Denotes  complex  conjugate 


Introduction 

Modern  aircraft  simulators  must 
represent  highly  nonlinear  vehicle 
aerodynamics  to  provide  effective 
pilot  training  and  accurate 
engineering  analysis.  However,  it  is 
still  useful  to  obtain  linear 
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aerodynamic  models  about  a  reference 
flight  condition  to  take  advantage 
of  powerful  linear  theory 
techniques,  especially  for  control 
law  analysis  and  specification 
compliance . 

The  Naval  Air  Warfare  Center's 
Manned  Flight  Simulator  (MFS)  has 
developed  a  highly  modular 
simulation  architecture.  The 
Controls  Analysis  and  Simulation 
Test  Loop  Environment  (CASTLE)  is  a 
generic  shell  structure  designed  for 
simulation  development,  execution, 
and  analysis. 

One  of  the  CASTLE  facilities  for 
engineering  analysis  is  the  MANeuver 
GENerator  (MANGEN) .  The  MANGEN 
facility  can  be  used  with  an  input 
data  file  to  overdrive  any  desired 
simulation  control.  Another  of  the 
CASTLE  facilities  is  the  Liner  Model 
Extraction  (LME)  program.  LME  is 
used  to  generate  a  linear  state- 
space  model  using  the  offset 
derivative  method. 

The  V-22  tilt-rotor  simulation  under 
the  CASTLE  architecture  was  run  in 
non-realtime  mode  to  produce  the 
simulation  data  for  the  linear  model 
comparisons.  MANGEN  was  used  to 
overdrive  the  control  input  for  the 
frequency  response  method,  and  LME 
was  used  to  extract  the  linear  model 
for  the  perturbation  method. 


7-22 _ simulation _ Overview 

The  V-22  Osprey  is  a  tilt-rotor 
aircraft  capable  of  flight  from 
hover  to  high  speed  airplane  mode. 
Control  of  the  V-22  is  accomplished 
through  both  conventional  airplane 
and  helicopter  controls.  The 
airplane  controls  include  elevator, 
rudders  and  flaperons,  and  the 
helicopter  controls  include 
longitudinal  cyclic,  lateral  cyclic, 
and  collective  pitch. 

The  V-22  has  a  digital  fly-by-wire 
control  system  consisting  of  a 
primary  flight  control  system  (PFCS) 
and  an  automatic  flight  control 


system  (AFCS)  .  The  PFCS  provides 
pilot  shaping  and  essential  feedback 
loops  for  primary  control.  The  AFCS 
is  designed  to  provide  level  1 
handling  qualities  and  auto-pilot 
functions . 

The  V-22  simulation  at  the  MFS  is  a 
high  fidelity  full  envelope 
simulation  utilized  for  pilot 
training  and  engineering  analysis. 
The  airframe  model  is  based  on  Bell 
Helicopter's  Generic  Tilt-Rotor 
Simulation. 


Initial  Conditions 

The  simulation  initial  flight 
conditions  were;  level  steady  state 
flight,  nacelles  fixed  in  airplane 
mode,  true  velocity  of  200  knots, 
and  constant  flap  setting.  The 
rotor  speed  and  flapping  controllers 
were  disabled  in  the  PFCS  to 
maintain  constant  swashplate 
positions.  The  core  AFCS  functions 
were  active.  Longitudinally,  this 
provides  a  pitch  attitude  command 
system  with  pitch  attitude  and  pitch 
rate  feedback  in  the  longitudinal 
stick  to  elevator  path. 


Perturbation _ Method 

The  Linear  Model-  Extraction  (LME) 
facility  under  the  CASTLE 
architecture  is  capable  of  obtaining 
linear  models  from  the  non-linear 
simulation,  in  the  standard  state- 
space  form: 

• 

x  =  Ax  -t-  Bu 
y  =  Cx  +  Du 

where 
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The  user  supplies  LME  with  the 
states,  state  derivatives,  inputs 
and  outputs  of  the  model  to  be 
extracted,  as  well  as  perturbation 
sizes  of  the  states  and  inputs.  LME 
uses  the  offset  derivative  method 
for  extracting  the  linear  models. 
In  this  method,  biases  are  added  to 
each  input  and  state,  and  the 
resulting  changes  in  outputs  and 
state  derivatives  are  recorded. 
Integration  is  frozen  within  the 
simulation  during  LME  so  that  the 
state  derivatives  are  not  allowed  to 
propagate.  A  schematic  of  this 
process  is  shown  in  figure  1. 
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The  state-space  matrices  are  then 


computed 
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Since  the  simulation  is  nonlinear, 
the  change  in  outputs  and  state 
derivatives  for  different  size 
perturbations  will  not  be  perfectly 
linear.  The  user  specifies  four 
perturbation  sizes  for  each  inp”.t 
and  state.  LME  perturbs  the 
simulation  positively  and  negatively 
for  each  of  the  perturbation  sizes 
and  picks  the  perturbation  size  that 
results  in  the  best  linearity.  A 
typical  example  of  the  perturbation 
selection  is  shown  in  figure  2. 


Figure  2.  LME  Perturbation 
Selection 


The  output  of  LME  is  the  state-space 
model  as  well  as  information 
describing  the  linearity  of  each 
coefficient . 

Model  Definition 

The  input  to  the  model  was  elevator 
position  and  the  outputs  were  pitch 
rate  and  normal  acceleration.  The 
states  chosen  to  represent  the 
aircraft  dynaunics  were: 
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Phase  (deg)  Mag  (dB) 


two  outputs  and  seventeen  states. 


X  body  velocity 
y  body  velocity 
z  body  velocity 
roll  rate 
pitch  rate 
yaw  rate 
roll  attitude 
pitch  attitude 
heading  angle 
altitude 


Since  the  rotors  contribute 
significantly  to  the  total  dynamics, 
the  following  states  were  included. 


left  longitudinal  flapping  angle 
right  longitudinal  flapping  angle 
left  lateral  flapping  angle 
right  lateral  flapping  angle 
rotor  speed 

right  gas  generator  speed 
left  gas  generator  speed 

The  resulting  model  has  one  input. 


Model  Reduction 

High  order  linear  models  are 
necessary  for  accurate  analysis  of 
control  systems.  Low  order  models 
are  often  desired  for  design  of 
control  systems  and  analysis  of 
flying  qualities.  In  this  case,  a 
fourth  order  model  representing  only 
the  longitudinal  dynamics  was 
desired. 

Some  of  the  states  in  the  extracted 
model  represent  off-axis  responses 
and  contribute  little,  if  any,  to 
the  longitudinal  dynamics.  The 
states  corresponding  to  the  rotor 
dynamics,  however,  affect  both  the 
phugoid  and  short  period.  These 
effects  need  to  be  included  in  the 
reduced  order  model. 
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The  model  reduction  method  chosen 
for  this  case  was  the  modal 
truncation  method.  Modal  truncation 
involves  transforming  the  model  into 
modal  form: 


X 


X  +  B  u 
m 


y 


c 


m 


X  +  D  u 
m 


where  A  is  a  mode  of  the  system. 
For  a  first  order  mode,  A  is  equal 
to  the  eigenvalue  (X)  .  A  second 
order  mode  is  represented  as: 

real(X)  imag(X) 

-imag(X)  real(X) 


The  contribution  of  each  mode  to  the 
overall  frequency  response  is 
calculated.  The  modes  with  the 
highest  contribution  are  retained, 
the  others  are  truncated. 


To  extract  a  fourth  order 
longitudinal  model  from  the  V-22 
simulation,  the  longitudinal  stick 
position  was  used  as  the  input 
signal,  and  the  desired  frequency 
responses  were  the  ub,  wb,  a^,  a^ , 
and  qb  responses  to  the  elevator 
position.  By  using  the  longitudinal 
stick  as  the  input  (instead  of 
directly  driving  the  elevator 
position)  the  control  system  reduced 
excursions  from  the  reference  flight 
condition.  Figure  4  is  a  simple 
schematic  of  the  simulation  process. 


ub,  wb. 


Figure  4  Schematic  of  Simulation 
Process 


The  phugoid  and  short  period  were 
easily  identified  with  the  model  in 
modal  format .  The  four  states 
corresponding  to  the  phugoid  and 
short  period  second  order  modes  make 
up  the  reduced  order  model .  The 
remaining  eleven  states  were 
truncated.  Figure  3  shows  the 
results  of  the  reduction. 


Frequency  Regponse  Method 

Overview 

The  frequency  response  method  uses  a 
time  history  input  which  contains 
the  frequencies  of  interest  to  drive 
the  simulation.  The  frequency  input 
and  the  output  of  the  simulation  are 
transformed  to  the  frequency  domain, 
and  the  desired  frequency  responses 
of  the  simulation  are  calculated.  A 
linear  model  is  then  fit  to  the 
frequency  responses  using  system 
identification  techniques. 


Input  Signal  Selection 

A  variation  of  the  swept  sine  input 
has  been  the  dominant  input  design 
for  frequency  response 
calculat ions . ^  ^ ^  The  design  goals 
for  the  frequency  input  were  to 
provide  good  spectral  excitement  in 
the  frequency  range  of  interest  and 
to  prevent  the  aircraft  (or 
simulation)  from  straying  too  far 
away  from  the  reference  flight 
condition . 

For  non-realtime  simulation  without 
a  pilot  in  the  loop,  a  swept  sine 
longitudinal  stick  input  drove  the 
V-22  simulation  away  from  the 
reference  flight  condition,  even 
with  the  AFCS  engaged.  Therefore,  a 
new  signal  that  prevented  large 
excursions  from  the  reference  flight 
condition  was  designed  using  a  sum 
of  sine  waves.  To  construct  the  new 
signal,  alternating  positive  and 
negative  sine  waves  at  selected 
frequencies  were  summed  into  a 
single  signal. 
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Frequency  (rad/sec)  Frequency  (rad/sec) 


(a)  Input  1  (b)  Input  2 

Figure  5.  Sum  of  Sine  Waves  Longitudinal  Stick  Inputs 


The  sum  of  sines  input  proved 
superior  to  the  swept  sine  input  for 
two  reasons.  The  low  frequency  input 
was  applied  throughout  the  entire 
signal  (instead  of  only  at  the 
beginning)  providing  better  low 
frequency  spectral  content,  and  the 
signal  was  easily  designed  to  limit 
the  excursions  from  the  reference 
flight  condition. 


Two  design  parameters,  the  frequency 
resolution  and  the  signal  length, 
were  considered  when  constructing 
the  sum  of  sines  input.  To  produce  a 
signal  that  has  equal  power  about 
the  half-period,  the  signal  length 
was  determined  from  the  input 
frequency  resolution. 


27I/A(i) 

^  t  _ 3  aw 

asw  ~  At  ~  At 

A  frequency  resolution  of  0.02 
radians  per  second  was  selected  to 
provide  good  low  frequency 
excitation  and  to  limit  the 
simulation  run  time.  Limiting  the 
simulation  run  time  reduced  the 
number  of  data  points  to  process  for 
the  frequency  response  calculation. 
The  length  of  the  sum  of  sines 
signal  was  calculated  using  a 
simulation  time  step  of  0.0333333 
seconds . 
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To  excite  the  V-22  fourth  order 
longitudinal  dynamics,  a  frequency 
range  from  0.05  to  10  radians  per 
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second  was  selected  based  on 
previous  simulation  data.  The  low 
frequency  of  0.05  radians  per  second 
captured  the  phugoid  dynamics,  and 
the  high  frequency  of  10  radians  per 
second  bounded  the  short  period 
mode . 

Two  input  signals  were  designed  to 
produce  data  from  two  separate 
simulation  runs.  Input  signal  1 
contained  45  sine  waves  with 
frequencies  selected  from  a  'vector 
of  0.06  to  10  radians  per  second. 
Input  signal  2  contained  45  sine 
waves  with  frequencies  selected  from 
a  vector  of  0.05  to  9.9  radians  per 
second.  The  frequency  vectors  for 
inputs  1  and  2  had  a  frequency 
resolution  of  0.02  radians  per 
second . 

The  frequencies  used  in  the  input 
signals  were  selected  to  be 
approximately  evenly  spaced  on  a 
logarithmic  scale  over  the  frequency 
range.  However,  not  quite  enough  low 
frequencies  were  available  in  the 
frequency  vector,  and  the  resulting 
signals  were  slightly  more  heavily 
concentrated  at  higher  frequencies. 

After  the  sine  waves  were  summed, 
the  signals  were  zero  padded  at  the 
beginning  and  end  to  create  a  total 
length  of  10,000  points.  The  zero 
padding  added  reference  flight  data 
to  the  beginning  and  end  of  the 
simulation  output. 

Finally,  the  signals  were 
appropriately  scaled  and  biased  for 
the  reference  longitudinal  stick 
position.  The  resulting  signals  are 
shown  in  Figure  5.  Notice  from  the 
power  spectral  densities  the  even 
power  content  at  the  discrete 
frequencies  contained  in  the  inputs. 

Simulation  Run 

The  V-22  simulation  was  initialized 
to  the  level  flight,  steady  state 
reference  condition,  and  the 
longitudinal  stick  inputs  shown  in 
Figure  5  were  used  to  drive  the 
simulation.  Overdriving  the 
longitudinal  stick  position  with  the 


input  data  files  was  accomplished 
using  the  CASTLE  MANGEN  facility. 

The  altitude,  pitch  attitude  and  u 
body  velocity  time  histories  for 
simulation  runs  1  and  2  are  shown  in 
Figure  6.  Notice  that  the  time 
histories  end  very  near  the 
reference  flight  condition  even 
after  a  simulation  run  of  more  than 
five  minutes. 

Frequency  Response  Calculation 

The  process  to  calculate  a  frequency 
response  from  frequency  sweep  time 
history  data  is  well  def  ined.  ^ ^ 
The  process  is  the  same  for  a  sum  of 
sines  input,  except  that  only 
discrete  frequencies  are  contained 
in  the  signal.  Therefore,  the 
simulation  output  contained  spectral 
information  only  at  the  discrete 
frequencies  included  in  the  input.  A 
frequency  bookkeeping  was  maintained 
when  processing  the  data  to  insure 
that  only  the  discrete  frequencies 
contained  in  the  input  were  used  to 
calculate  the  resulting  frequency 
responses . 

The  ub,  wb,  a^,  az<  time 

histories  from  simulation  runs  1  and 
2  were  concatenated,  and  the  data 
was  detrended  to  remove  the  bias  and 
linear  drift.  The  concatenated  data 
was  windowed  with  two  non¬ 
overlapping  10,000  point  rectangular 
windows.  Large  window  sizes  give 
better  low  frequency  resolution,  and 
small  window  sizes  average  the  data 
to  reduce  high  frequency  noise. 
Since  the  simulation  time  history 
data  contained  essentially  no 
noise,  minimal  averaging  was 
necessary . 

The  chirp  z-transform  was'  used  to 
process  the  time  history  data.  This 
advanced  Fourier  transform  provided 
a  much  finer  frequency  resolution 
than  the  standard  fast  Fourier 
transform.  The  frequency  resolution 
of  the  chirp  z-transform  is 
determined  from  the  following 
equation. 
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UB  (ft/sec)  Pitch  Attitude  (deg)  Altitude  (ft) 
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(10.0  -  0.0) 
10000 


0.001 


The  chirp-z  transform  function 
included  in  the  MATLAB  Signal 
Processing  Toolbox  was  used  to 
transform  the  time  history  data  to 
the  frequency  domain. 


indicates  than  none  of  the  output  is 
caused  by  the  input .  The  coherence 
is  calculated  from  the  following 
equation . 
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Once  the  data  is  transformed  to  the 
frequency  domain,  the  input  power 
spectral  density,  the  output'  power 
spectral  density,  and  the  input- 
output  cross  spectral  density  were 
estimated  from: 
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The  single  input  single  output 
transfer  function  was  then 
calculated  from: 


H  (M) 
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After  the  transfer  function  was 
calculated,  the  magnitude  in 
decibels  was  found  from: 


The  frequency  bookkeeping  maintained 
when  calculating  the  frequency 
responses  for  the  concatenated  data 
produced  magnitude,  phase  and 
coherence  data  at  90  discrete 
frequencies.  The  data  were  slightly 
more  heavily  concentrated  at  higher 
frequencies  than  at  lower 
frequencies.  To  produce  data 
approximately  evenly  spaced  on  a 
logarithmic  scale,  52  of  the  90 
discrete  frequencies  were  selected. 


State-space _ Model _ _ L2 _ Ul£ 

Frequency  Response 

Once  the  frequency  responses  were 
calculated  from  the  simulation  data, 
a  fourth  order  longitudinal  state- 
space  model  was  fit  to  the  data.  The 
four  states  of  the  model  were  ub, 
wb,  qb,  and  0;  the  single  input  was 
5^:  and  the  outputs  were  ub,  wb,  ax, 
a^f  and  qb. 

The  following  equation  is  used  to 
calculate  the  frequency  response  of 
a  state-space  model. 

H(a))  =  C[(jQ))I  -  a]"1b  +  D 


mag  =  201ogl0  (  |  H  (OJ)  |  ) 
s  sw 

and  the  phase  in  degrees  was  found 
from: 

180 

phase  =  - imag  ( log  (H  (to)  )  ) 

ssw  JIC 

The  coherence  function  is  a  measure 
of  how  much  of  the  output  is 
linearly  related  to  the  input.  A 
coherence  of  one  indicates  that  all 
of  the  output  is  produced  by  the 
input,  and  a  coherence  of  zero 


A  cost  function  was  set  up  to 
calculate  the  magnitude  and  phase 
error  between  the  fourth  order 
state-space  model  and  the  sum  of 
sine  waves  frequency  responses.  The 
cost  function  was  designed  so  that 
seven  decibels  of  magnitude  produced 
the  same  error  as  one  degree  o* 
phase . 
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Solid  =  Identified  State-Space  Model  o  =  Frequency  Response  Points 


Frequency  (rad/sec) 


Flgur*  7.  Identified  Model  Fit  to  Frequency  Responses 


52 

cost  =  X  (7  (mag  {(1)^)  -  mag  (o)^ )  ^^)  )  2 

i  =  l 
52 

+  S  (phase(O)^)  -  phase  (oa^  )  _  2 

i  =  i 

The  index  i  is  used  to  sum  the  error 
at  each  of  the  52  discrete 
frequencies.  The  cost  function  was 
minimized  using  a  binary  search 
method  called  digital  matching.^ 


The  state-space  model  fits  to  the 
pitch  rate  to  elevator  and  normal 
acceleration  to  elevator  frequency 
responses  are  shown  in  Figure  7  . 
Every  other  frequency  point  used  in 
the  identification  was  plotted  for 
clarity.  Notice  the  coherence  for 
the  frequency  responses  are  very 
near  one  throughout  the  frequency 
range  of  interest. 
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Solid  =  Non-Linear  Simulation  Dashed  =  Linear  Models 


Plffur*  8.  Linear  Model  Time  History  Verification 


- VTlflcatlon 

The  state-space  models  extracted 
using  the  perturbation  and  frequency 
response  methods  were  both  verified 
against  a  time  history  response  from 
the  non-linear  simulation.  The 
nonlinear  simulation  was  driven  with 
an  elevator  doublet  about  the 
reference  flight  condition.  Both 
linear  models  were  then  driven  with 
the  same  elevator  doublet.  The 
linear  model  responses  compared  to 
the  simulation  response  are  shown  in 
Figure  8. 


Modal _ Comparison 

The  short  period  and  phugoid  modes 
of  each  of  the  models  were  compared 
by  calculating  the  undamped  natural 
frequency  and  damping  ratio  of  the 
two  modes  from  the  eigenvalues  of 
the  system  matrices.  Table  1 
contains  the  undamped  natural 
frequency  and  damping  ratios  for 
each  method . 

The  frequency  response  was 
calculated  from  the  perturbation 
method  linear  model  and  compared  to 
the  frequency  response  from  the 
frequency  method.  The  results  are 
shown  in  Figure  9 . 
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Solid  =  LME  Model  Dashed  =  Frequency  Response  Method 


Flgurs  9.  Frequency  Response  Comparison  of  Perturbation 
and  Frequency  Response  Methods 


Short 

Period 

Phugoid 

COn 

(On 

perturb. 

method 

2.01 

0.497 

0.165 

0.192 

frequency 

response 

method 

2.07 

0.513 

0.168 

0.187 

Table  1.  Undamped  natural  frequency 
and  damping  ratio  for  the  short 
period  and  phugoid  modes. 


Concluaions 

Fourth  order  longitudinal  state- 
space  models  were  extracted  from  the 
V-22  full  nonlinear  simulation  using 
the  frequency  response  and  the 
perturbation  methods.  Time  history 
verification  against  the  full 
nonlinear  simulation  showed  that 


both  models  accurately  predicted  the 
pitch  rate  and  normal  acceleration 
response  to  an  elevator  doublet . 
Comparison  of  the  short  period  and 
phugoid  modes,  as  well  as  the 
frequency  response  of  the  models, 
showed  that  similar  dynamic 
characteristics  were  identified  by 
both  methods .  Within  the  scope  of 
this  paper,  both  the  frequency 
response  and  perturbation  methods 
were  effective  means  to  extract  a 
low  order  state-space  model  from  a 
nonlinear  simulation. 

R^eomm^ndat Iona 

The  primary  advantage  to  the 
perturbation  method  is  that  the 
linear  model  is  computed  directly. 
The  frequency  method  requires 
identifying  a  model  to  fit  the 
frequency  data,  which  can  become 
significantly  more  difficult  with 
higher  order  models.  When  higher 
order  models  are  desired,  the 
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perturbation  method  offers  the  most 
direct  solution. 

In  some  cases,  as  in  classical 
stability  analysis  of  feedback 
loops,  the  frequency  response  is  the 
desired  end  result.  The 
perturbation  method  requires 
computing  the  frequency  response 
from  the  linear  model.  Although 
this  was  a  simple  task  for  a  fourth 
order  model,  computing  the 
frequency  response  for  very  high 
order  systems  can  become 
computationally  difficult.  The 
frequency  method  computes  the 
frequency  response  directly  and  is 
unaffected  by  the  order  of  the 
system.  When  the  frequency  response 
of  large  systems  is  desired,  the 
frequency  method  should  provide  the 
most  efficient  solution. 
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